/*
 * Bisection.cpp
 *
 *  Created on: Aug 26, 2010
 *      Author: snirgaz
 */

#include "Bisection.h"


void Bisection::init(WorldLines *wl, QmcParameters * qmcParams, MoveParams *moveParams,
		MTGen *RndGen) {
	wl_ = wl;
	qmcParams_ = qmcParams;
	moveParams_ = moveParams;
	int M_ = lat_->wlData_.M_;
	zeroOne_.init(RndGen, UniformReal::param_type(0., 1.));
}

bool Bisection::generateBBMove() {
	energyDiff_.zero();
	if (this->generateBB())
		return true;
	else
		return false;
}
bool Bisection::generateBB() {
	bool flag;
	// Iterate Levels
	flag = true;
	while ((moveParams_->beads.size() < moveParams_->mbar) && (flag))
		flag = this->nextLevel();
	return flag;
}
bool Bisection::generateRWDMove() {
	int numOfBeads;
	double p;
	bool flag;
	EnergyData ED;
	moveParams_->currentBeads.clear();
	// Generate Jump
	numOfBeads = moveParams_->mbar;
	BeadVec::const_iterator beadStart = moveParams_->beads.begin();
	BeadVec::const_iterator beadEnd = --moveParams_->beads.end();
	BeadVec::iterator beadTau = boost::next(moveParams_->beads.begin(),
			moveParams_->tau);
	// Add Jump Bead
	randomWalks_.generateJumpBeadPos(beadStart, beadEnd, beadTau);
	// Check If Bead Tau is a jump
	beadTau->proposedJump_ = true;
	if ((beadTau->WLSlice_ == wl_->wlData_.wormSlice_)
			&& (wl_->wlData_.isOpenClosed_ == OPEN))
		beadTau->originalJump_ = true;
	moveParams_->currentBeads.push_back(moveParams_->tau);
	assert(
			(beadTau->RWSlice_ > 0)
					&& (beadTau->RWSlice_ < (moveParams_->mbar - 1)));
	// Energy Calc of Jump Bead
	ED = energyCalc_->calcBeadsEnergy();
	p = exp(-ED.getEnergy());
	assert(isInfinite(ED.energy_));
	if (!(accept(p)))
		return false;
	moveParams_->energyDiff.update(ED);
	// Left and right parts
	// Left
	flag = this->generateBB();
	return flag;
}
void Bisection::nextLevelBeadPos() {
	list<int>::iterator itb, itb_next;
	int deltaPos, i, numOfNewBeads;
	double ratio;
	ParticlePos meanNewPos, start, end, delta;
	// Clear Current beads Data
	moveParams_->currentBeads.clear();
	moveParams_->deltas.clear();
	itb = moveParams_->allBeads.begin();
	// Generate New Locations
	while (boost::next(itb) != moveParams_->allBeads.end()) {
		// FirstBed
		Bead &aBead = moveParams_->beads[*itb];
		// Second Bead
		Bead &bBead = moveParams_->beads[*boost::next(itb)];
		deltaPos = aBead.RWSlice_ - bBead.RWSlice_;
		// Check if adjucent beads
		if (deltaPos != 1) {
			// New Bead
			Bead &beadCurrent =
					moveParams_->beads[aBead.RWSlice_ + deltaPos / 2];
			beadCurrent.RWSlice_ = aBead.RWSlice_ + deltaPos / 2;
			beadCurrent.WLSlice_ = (aBead.WLSlice_ + deltaPos / 2) % M_;
			beadCurrent.n_ =
					(beadCurrent.WLSlice_ > moveParams_->mStart) ?
							(moveParams_->wlStart) : (moveParams_->wlEnd);
			if ((lat_->wlData_.isOpenClosed_ == OPEN)
					&& (beadCurrent.WLSlice_ == lat_->wlData_.wormSlice_))
				beadCurrent.originalJump_ = true;
			else
				beadCurrent.originalJump_ = false;
			// Set Bead Mean pos
			if (beadCurrent.proposedJump_)
				start = beadCurrent.jumpPos_;
			else
				start = beadCurrent.beadPos_;
			end = bBead.beadPos_;
			delta = end - start;
			//end.subBeadPos(delta, start);
			ratio = double(deltaPos / 2) / double(deltaPos);
			delta = ratio*delta;
			//delta.mulBeadPos(ratio);
			meanNewPos = start + delta;
			//start.addBeadPos(meanNewPos, delta);
			beadCurrent.beadPos_ = meanNewPos;
			// Delta List
			moveParams_->deltas.push_back(deltaPos);
		}
		itb++;
	}

}
void Bisection::nextLevelRW() {
	randomWalks_.generateLevelPath(moveParams_);
}
bool Bisection::nextLevel() {
	double p;
	EnergyData energy;
	this->nextLevelBeadPos();
	this->nextLevelRW();
	energy = energyCalc_->calcBeadsEnergy();
	p = exp(-energy.getEnergy());
	assert(isInfinite(energy.energy_));
	if (accept(p)) {
		moveParams_->energyDiff.update(energy);
		return true;
	} else
		return false;
}

bool Bisection::calcSwap(int originalWL, int swapToWL) {
	EnergyData energySwap = energyCalc_->calcEnergySwap(originalWL, swapToWL);
	double p = exp(-energySwap.getEnergy());
	assert(isInfinite(energySwap.energy_));
	if (accept(p)) {
		moveParams_->energyDiff.update(energySwap);
		return true;
	} else
		return false;
}
